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A higher order panel method has been used to simulate 
three dimensional wake roll-up for large aspect ratio symmetric 
wings in steady inviscid incompressible flow. 

Two methods have been considered for this purpose. The 
first method v;hich will be referred to as vorticity distribution 
method - consists of determination of wake geometry iteratively 
in three steps. The first step aims at computation of linear 
vorticity distribution over lifting surface using appropriate 
boundary conditions. The second step is concerned with determi- 
nation of induced velocities. This is followed by final step 
of updating wake geometry. These steps are performed repeatedly 
to achieve convergence. 

This method leads to certain numerical difficulties which 
have been alleviated in second method here after called 
circulation distribution method by replacing linear vorticity 

i 

by quadratic circulation distribution over the lifting surface, 
other steps remaining unaltered. This has been successfully 
employed for generation of wake roll-up geometry for rectangular 
wings in steady flow, : 

The vorticity distribution method could only be used for ' 

'i 

obtaining the vorticity distribution on the wing surface and * 

hence the pressure distributions over a flat plate rectangular ; 

•I,. 

wing of large aspect ratio. The circulation distribution method 
was used to obtain the wake roll up on large aspect ratio flat \ 
plate wings at angle of attacks upto 10®, J 

I 

. , .... , . ..... .. ..... ....... ... . ...... . .. ... . .... .2 





CHAPTER - 1 


IRTRODUCTIOn 


1.1 Preliminary remarks 

The knowledge of the trailing vortex wake is required 
because of the increasing importance of interactions beta^reen 
aircraft in crowded air corridors and the need to minimize 
these interactions by accelerating the dissipation of the 
trailing vortices. 

A wing of finite span in flight generates trailing 
vortices which rollup into two "discrete" vortices. During 
the rolling up process the trailing vorticity is displaced 
vertically downwards, assuming that the aircraft is in 
horizontal flight, under the influence of the wing circulation 
and the self-induced downwash velocities, 

A starting point for the calculation of the trailing 
vortex wake is the realistic evaluation of the vorticity 
distribution on the wing and in the near wake. The calculation 
of a vorticity distribution on the wing planform and its 
trailing wake leads to an Important result? the accurate 
evaluation of wing aerodynamic characteristics including 
the nonlinear effects. 

Fob the vortex flow behind a wing of finite span at 
angle of attack, there exists a pressure differential 
tween the top and bottom surfaces as a bi^roduct of the 
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lift force. The vorticity is shed as a result of the 
difference in the direction of the flow between the wing 
surfaces at the trailing edge. 

For high Reynolds number flow past a thin wing at 
moderate angles of attack, convection dominates diffusion 
and the vorticity is confined to a thin, free shear layer; 
under the influence of the self induced velocity, the shear 
layer has tendency to roll-up into vortex cores. In the 
downstream direction the vorticity is continuously fed into 
vortex cores resulting in a growth both in strength & dimension. 
Further downstream of the wing the roll up process will be 
completed and most of the vorticity is contained within the 
cores. At this stage the velocity gradients in the cores are 
large and the viscous effects beccxne more important, conseguently 
instability may occur which can result in the break up of the 
trailing vortex system. We are considering only initial stage 
of the wake roll-up - potential flow is considered with embedded 
vortex sheets. Thus we are able to predict viscous phenomenon 
using potential solution without solving Navier-strokes equations 

1,2 literature Review 

In this section earlier attempts to model lifting-surface 
flows are reviewed and evaluation of the present method is 
discussed. 
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The earliest analytical method was presented by Prandtl 
in 1918 (See e«g, Karamcheti /2/) and is known as the lifting , 
line theory. He considered rectangular wings having large 
aspect ratio at small angles of attack. The bound vortex sheet 
was approximated by a single bound vortex line having variable 
strength. Spatial conservation of vorticity requires a 
continuous sheet of trailing vortices to be shed from the bound 
vortex line. The resulting model in assumed to lie in a plane 
parallel to the free stream. In addition the model was 
simplified by assimptions concerning the spanwise component 
of the bound vortex line. The model is very simple and limited 
to rectangular wings having large aspect ratios and low angles 
of attack. 

Using the method of matched asymptotic expansion. Van Dyke 
/14/ and Thurber /I 6/ refined Prandtl ‘s analysis to obtain 
better approximations, Kemey /36/ corrected Van Dyke's expressior 
for the lift-curve slope. However these expansions possess a 
number of non-uniformities near the wing tips, discontinuities 
in planform, and root juncture of a swept back wing among others. 

One of the earliest models for rectangular wings having 
low aspect ratios was presented by Bollay /5/, In the limiting 
case of a rectangular wing of zero aspect ratio (finite span 
but infinite chord) , he assumed that the wing surface is 
equivalent to a distribution of spanwise, bound vortices, ®ie 
trailing vortices are shed rearwards from the wingtips at half 
the angle of attack in a plane perpendicular to the wing surface. 
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Assrair.r the lift distribution and the dovm\';ash to be constant 
across the span, he expressed the total normal force coefficients 
as a nonlinear function of the angle of attack. He also extended 
the theory to v/ings having finite but small aspect ratios. Here, 
the trailing vortices from the vring tips are again straight taut 
follow the direction of the resultant velocities at the wing 
tips. The no-penetration boundary condition v/as satisfied along 
the mid-chord, leading to an integeral equation. The form of the 
function giving the vorticity in terms of the distance frx>m the 
leading edge v/as assumed. This method yields surprisingly good 
results for wings having aspect ratio less than one, hat. it does 
not predict the distribution or shape of the trailing vortices, 
and it cannot be applied to other planforms, 

Multhopp /&/ presented a method, based on the linear potential! 
equation obtained by the use of continuity and momentum equations I 

i 

for steady inviscid flow. This method uses fields of doublets over 1 

i 

the projection of the wing on a plane parallel to the free stream. * 

f 

The downwash velocity condition is satisfied on the projection * 

of the wing. This results in an integral equation relating the i 

local-incidence and the distributed-load coefficient. The Kutta \ 

I 

condition is satisfied by assuming that the load functions [ 

disappear toward the trailing edge. Finally the integral equation | 
is solved by selecting a limited number of independent load : 

distributions and satisfying the integral equation by a linear ,, 

I 

combination of these distributions at a certain number of so i 

t 

f 

called '"pivotal points". This method is limited to small angle 
of attack and large aspect ratio. 
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Similar methods were developed to model flow over very 
highly svrept wings. These were based on the slender-body 
approximations. One such method is due to Lawrrence /!/» who 
also presented an extensive summary of previous methods. 

Lawrence did not attempt to model separation at sharp leading 
edges and/or wing tips, which can occur even at moderate angles 
of attack. 

Apparently, Legendre made the first attempt to include 
leading-edge separation in his model of flows over delta wings. 

In his model, isolated vortices were placed over the wing. Their 
strength and position were determined by forcing the flow to be 
tangential to the wing surface at the leading edge. Brown and 
Michael /8/ improved this model by adding feeding vortex sheet 
extending from the leading edge to the isolated vortices above 
the wing. The strengths and positions of the vortices are obtained 
by requiring the flow to be tangential at the leading edge and 
the resultant force on the feeding sheets and the vortices to 
be zero. They did not make the pressure difference across the 
feeding sheets zero everywhere. Mangier and Smith /9/ modified 
this model further by assuming a curved shape for the feeding 
vortex sheets, and by forcing the pressure difference across 
these sheets to vanish at selected points. In these methods the 
shape of the vortex sheets at the leading edge is prescribed, 
the flow is assumed to be conical, and the slender body assxamption 
is used. Consequently these methods are restricted to slender 



delta wings and most importantly small angle of attack. These 
assxamptions do not appear to be justified for the delta wings 
used in the wind tunnel tests. 

In an attempt to treat wings having general planforms 
and arbitrary aspect-ratios, Gersten asstiraed that trailing 
vortices are shed from each point of the lifting surface at one 
half of the angle of attack. Such an approach leads to a non- 
linear relationship between the total loads and the angle of 
attack. Gamer and Lehrian /13/ combined the method of Multhopp 
and Gersten, Their model is based on essentially the same assTomptio 
which were used by Gersten, 

Theoretical and semi-empirical methods were presented by 
Sacks et al /23/, The methods are based on the assiimption that 
the lift coefficient is the stun of two terms. The first is the 
lift coefficient obtained by the method of Lawrence, vdiile the 
second is due to the contribution of the leading edge separation. 
They considered delta wings of small aspect ratios and hence 
used the slender body approximation. The vortex sheet shed from 
the leading edge is represented by a finite number of vortex 
filament. Here, the basic unknowns are the strengths and position 
of these vortex filaments. One equation was obtained by satis- 
fying the Kutta condition at the leading edge 'vAile the other was 
obtained from shedding rate. The shedding rate was determined 
by either a semi-CTiperical or a theoretical method. The results 
of theoretical method over predicted the lift coefficient, 
while those of the s®ni-«npirical method gave satisfactory results, 
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Pinkleman /35/ extended this method and obtained a solution 
for the interacting Canard-wing system. Other methods were 
developed on the basis of treating the shape of the trailing 
vortices by the two dimensional analysis associated with the 
intersection of the vortex sheet and the Trefftz plane 
Rossow /44/. 

Polhamus /31/ developed a method to predict the total 
lift and the drag due to lift for sharp edge delta and delta- 
like planfoms. The method is restricted to thin wings having 
no camber and no twist. In addition the leading edge is assumed 
to have sufficient sharpness so that the separation is fixed 
at it and no leading edge suction is developed, 

Polhamus /30/ extended this method to account for 
compressibility by using the Prandtl-Glauert transformation. 
Bradley# et al /45/ extended Polhamus 's method to more general 
planfoxms. Although Polhamj^s’s method predicts the total lift 
coefficient with sufficient accuracy# it is incapable of 
predicting the pitching moment and the shape of the wake, 

Kuo Sc Morino /43/ presented a method for determining the 
distributed and total aerodynamic loads on lifting bodies 
having arbitrary shapes and sharp trailing edges in a steady# 
subsonic potential flow. The Prandtl Glauert transformation 
is used to reduce the linearized compressible flow equation 
to Laplace's equation. Green's theorem is used to obtain an 
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integral equation relating the velocity potential to its 
normal derivative on the surface enclosing the lifting body 
and its wake at the trailing edge. The shape of the wake is 
prescribed. Morino and Kuo /47/ extended this method to account 
for xinsteadiness in the flow. 

In brief, the most common drawback of all the analytical 
methods presented above is that none adequately models the 
wake. Hence they are restricted to small angles of attack. In 
addition they are restricted to wings having certain planforms, 
aspect ratios, and/or small thickness ratio. 

Discrete-vortex methods are based on approximating the 
continuous distribution of the bound vorticity (representing 
the lifting surface) by either a finite number of horseshoe 
vortices or a vortex lattice composed of a finite number of 
vortex filaments and the continuous distribution of free 
vorticity (representing the wake) by a finite nxamber of vortex 
lines which extend downstream. This technique has been known 
for a long time. The earliest models were presented by Prandtl 
/2/ and Bollay /5/ but their further development and refinement 
only became practical with the advent of digital ccatputers. Ihese 
methods can be divided into two types. Linear and non-linear 
discrete vortex methods. 

Numerical methods of the LINEAR DiSCRETS-VORTEX type were 
presented by several investigators such as Rubbert /1 2/, 

Hedtpan /20/ and Belotserkovskii /2l/, for solving the steady 
flow probl«n Later these methods were extended to the unsteady- 
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flow problem (oscillatory case) by combining the vortex-lattice 
with a doublet-lattice by Kalman# et al /22/# Albano and 
Rodden /27/, Rodden# et al /38/, Giesing et al /36/ Ashley and 
Rodden /34/# Hua /4l/ and Rodcien /49/« Only steady flow problem 
will be discussed here. 

In this approach# the so-called vortex lattice method 
(VLK) # the lifting surface is replaced by a number of panels 
arranged in strearawise columns. Each panel is replaced by a 
horseshoe vorteit with its bound part positioned on the quarter 
chord line of the panel, Tj^ich is the location of the aerodyna- 
mic centre of a two dimensional plate in a steady incompressible 
inviscid flow. With this arrangement a vortex lattice represents 
the lifting surface. The wake attached to the lifting surface is 
represented by rectilinear-vortex lines which extend downstream 
to infinity. The direction of these lines is prescribed to be 
parallel either to the free stream velocity or to the root chord 
of the lifting surface. With this configuration# the only unknown 
is the circulation distribution. This is obtained by satisfying 
the no penetration boxindas^^ condition at a certain number of 
points (known as control points) on the wing surface. The control 
points are centered spanwlse at the % chord length of each 
panel. The reason for .the ' selection of this location for control 
point is based on numerical experiments, James /28/ and Hough /46/ 
indicated that this location yields fast convergence of the 
aerodynamic loads as the number of panels is increased. Once the 
circulation distribution is obtained the forces due to the botmd 



vortices can be calculated. It is noted that a bound 
vortex line runs betv?een the leading row of control points 
and the leading edge of a delta iv'ing and betv/een the out board 
columns of control points and the wing tip for other wings. 

Thus separation at the leading edge and v/ing tip is not taken 
into account. Hence it is restricted to small angle of attack. 

It is worth mentioning that seme investigators introduced 
compressibility effects into this technique by using Prandtl- 
Glauert transformation e.g. Hedman /20/ and Kalman et al /22/, 

Hess /39/ used a finite-strength vorticity distribution 
over the wing surface instead of vortex lines. He used a step 
function for the spanwise variation of vorticity# while he 
used a constant strength around each chordwise strip. This 
method is also limited to small angles of attack. 

All the previous models suffer from a caranon basic fault# 
none of them models the leading-edge and wing-tip vortex system 
adequately. For the intended applications, these vortex systems 
strongly influence the flow in the immediate neighbourhood of 
a significant portion of lifting surface; thus methods whidh 
do not describe vortex systems accurately cannot be expected 
to predict the aerodynamics characteristics of the wing accuratel^j 

In the NON-LINEAR DISCRETE VORTEX type^the first attempt 
to model the wake was made by Srmolenko /1 9/, who considered 
rectangular wings having small aspect ratios. He replaced the 
wing surface by a system of spanxvlse vortex lines. The circulatior 
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are constant along each vortex line but vary from line to line. 
The vortices in the wing tip wake are considered to be recti- 
linear and to lie in two vertical planes which are perpendicular 
to the lifting surface. Each trailing vortex is directed parallel 
to the projection of the resultant velocity onto these planes 
at the point where the spanwise lines intersect the wing tips. 

The boundary conditions are satisfied on the average along 
spanwise lines distributed at the three-quarter chord lines 
of each panel. This method is subject to the same limitations 
as Bollay's method /5/. 

Belotserkovskii /25/ greatly improved Ermolenko’s method 
and developed the first successful model of wing- tip wake. He 
also replaced the wake by a finite ntxmber of discrete vortex 
lines / but in contrast with the earlier methods, he made each 
line consist of a finite number of segments, where the last 
one is semi- infinite segment. Then through an iterative 
procedure, he aligned each finite segment with the velocity 
at its mid-point. He showed the method to be highly reliable 
in predicting the total and distributed loads on parallelogram 
wings at large angle of attack. 

Butter and Hancock /33/ replaced the wing by a single 
lifting-line bound vortex which represent the circulation ' 
distribution around the wing while the trailing vortex sheet 
is approximated by a number of discrete line vortices. Starting 
with a planar system and following step-by-step process, working 
downstream aft of the wing trailing edge, the deformation of 
the trailing vortex sheet is determined. 
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Mook and Maddox /48/ used a vortex lattice method and 
considered leading separation. They modified the steady portion 
of the program of Giesing, Kalman and Rodden /36/ by super- 
posing a system of discrete vortex lines on the existing 
lattice. They determined the shape of the lines trailing from 
the leading edge. However the method cannot account for wing 
tip and trailing edge vortex sheets being free of forces. 

Through the advance of computer technology, a number of 
fully three dimensional models have been developed for the 
prediction of wake roll-up. The vortex-lattice method has foimd 
wide applications for canputing the attached flow loading of 
lifting surface. Rom and zorea /42/ developed a computational 
procedure to relax the wake. They use finite number of line 
vortices to represent the trailing vortex sheet. The wake 
geometry is evaluated in an iterative fashion starting from an 
assumed initial wake position. Sucio and Morino /51/ (developed 
a method that uses constant doublet panels on the wing and the 
wake to calculate the wake roll-up behind rectangular wings. The 
use of constant doublet panels for the wake modelling is equivalen 
to using vortex filaments. In recent years the constant doublet 
panel method has been widely applied to helicopter applications 
to model the blade and wake /57/, 

For the cases when the vortex wake passes close to lifting 
surface, the filament wake model may not be adequate to describe 
the flow details. This is true for leading edge sheets of low 
aspect ratio wings. It is also noted that a hi^er order panel 
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method solution for the wake roll-up behind a large aspect 
ratio wing is not available in literature, 

Yeh /61/ used linear vortex panel method to calculate 
the three dimensional wake roll-up behind large aspect ratio 
symmetric wing in steady inviscid incompressible flow. He 
obtained the wing spanwise circulation distribution from 
lifting line theory and wake geometry was evaluated in an 
iterative fashion. The capability of the method to model a 
rectangular wing with a deflected flap was also studied. 

Higher order panel methods for three dimensional vortex 
calculation have been developed by Johnson et al /54/, 

Hoeijmakers and Co-workers /59/ and Kandil et al /60/, These 
methods have been used to study leading and side edge separation 
for low aspect ratio wing. The method developed at Boeing /54/ anc 
NLR /58/ use second order doublet distribution. 

In the method of Kandil et al /60/, vortex panels 
(triangular panels in the wake) with a first order vorticity 
distribution are used in the near field calculations. In the 
far field calculations, the distributed vorticity is Itimped 
into equivalent concentrated vortex lines. We choose to modify 
this method and apply it to large aspect ratio wing. 

An excellent review of the literature on vortex wakes 
is given by Hoeijmakers /59/. Many of the studies use a two- 
dimensional time dependent analogy to study the wake roll-up 
for large aspect ratio case. The technique is d^eveloped under 
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the assumption that the flov? varies in the streamwise 
direction and that there is no direct influence due to the 
upstream v/ing. The two dimensional time dependent problem can 
be considered as a marching problem in the sense that vortex 
sheet is determined at each cross flow plane and is convected 
with its own induced velocity in subsequent cross flow planes. 

The work by Pink and Soh /53/ appears to be most successful 
in producing smooth wake roll-up by using the discrete vortex 
method. Murman and Stremel /56/ utilized the “Cloud- in Cell" 
method to capture the vortex wa]ce motion, A sophisticated 
second order doublet panel method (curved segments) developed 
at NLR /59/ is capable of describing canplicated vortex sheet 
motion and vortex cores are also included to model highly 
rolled-up regions. 

All users of the vortex representation have met difficulties 
of convergence i ,e, growing randomness of location of equivalent 
vortices so that the line Joining these crossed itself and could 
no longer serve to represent a sheet. This inconsistency appeared 
to become more pronounced as segmentation was refined. Brikhoff 
/ll/ asserted the inevitability of such randomization as a 
consequence of the Hamilton principle applied to a two-dimensional 
systemi of discrete vortex filaments. Von Karman /l/ had already 
d^nonstrated strongly growing instability with decreasing 
separation for the special case of linear array of point vortices 
of equal strength* 
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A suitable problem for testing a new approach is that 
of the calculation of roll-up of the trailing vortex sheet 
well behind an elliptically loaded lifting surface. This has 
challenged many v/riters. Westwater /4/ x^as the first to use the 
Roshenhead /3/ method for this problem. He confined his attention 
to points far enough downstream of the wing for the influence 
of the equivalent bound vortex to be negligible, thus converting 
the problem to a two dimensional one; namely a trefftz plane 
calculation. Fortuitously, Westwater used only ten discrete 
vortices, a small enough number to escape serious irregularity 
of vortex position in the limited period of his calculation. 
Attempt to reproduce Westwater* s result by Takami /15/ and 
Moore /32/ were not successful. For much finer subdivision than 
that used by Westwater, these authors found that the vortices 
moved Chaotically and that the sheet or rather the straight 
line joining appropriate vortices, crossed through Itself, 

Moore /32/ showed that this phenomenon was not due to round 
off errors. If at scane instant the vortex strength was reversed 
the chaotic motion would unscramble and the vortices would return 
to their original position. Several authors then put forward 
schemes intended to overcome the difficulty, Moore /50/ used 
a different scheme. If the vortices near the tip are too far 
or the distance between the vortices is greater than critical 
distance, they are combined to form a tip vortex or line vortex 
with its coordinates at the centroid of the tat^o vortices and the 
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circulation is equal to the combined circulation of original 
pair. We have used this method- for our problem, 

1.3 Present Scope 

In the present vrork an atterapt is made to calculate the 
wake roll-up behind large aspect ratio wings, for the accurate 
computation of aerodynamic coefficients over the lifting 
surfaces, 

# 

Vorticity distribution over the panels in the near wake 
is determined, and then induced velocities are calculated at 
the comer point or average point of triangular panels chosen 
to model the wake. These velocities have been used to obtain 
the shape of free wake. Now aerodynamic coefficients are 
calculated using these new coordinates for wake surface rather 
than a planar wake. 

The present wake model bears no small perturbation 
assumptions and all of the equations engaged are solved 
analytically. Moreover the method is capable of calculating 
the wake geometry for a variety of planforms and shapes when 
spanwise circulation distribution is known either from numerical 
calculation or experiment. 

The first chapter deals with the requirement of trailing 
vortices, scmae attenpts by various authors to model the lifting 
surface flows and a •brief summary of the probl«n considered is 
discussed. 

In the second chapter the axes sys-tem for -the problem, 
mathematical formulation or governing equations and boundary 
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conditions are discussed. 

The third chapter discusses the method of solutions for 
the wake rollup problem and the manner in which the boundary 
conditions are satisfied, 

Fourtli chapter describes the computer program detailing 
the function subprograms and subroutines, indicating the various 
inputs and outputs, A list of variables used in the program is 
also attached. 

The last chapter gives the results of canputation and 
discussion of these results. 



CHAPTER - 2 


PROBLEM SPECIFICATION AND I/LATHEMATICAL FCRIX'LAriCI-: 

The accurate calculation of the position and strength 
of the vortex wake behind a blade in motion is a critical 
problem for rotorcraft applications. A detailed model of the 
vortex wake would permit substantial improvements in many 
aspects of rotary wing design such as performance, acoustics, 
vibration and response. Trailing vortices are of interest to 
minimize the interaction between the two aircraft flying in 
vicinity with each other. Purpose of the present investigation 
is to canpute the wake roll-up behind large aspect ratio wing, 
for the accurate calculation of aerodynamic coefficients over 
the lifting surfaces. 

Consider the flow of a steady freestream of speed Uo«» at 
an angle of attack a past a thin, wing of large aspect ratio. 
Since the flow field considered here is the high Reynolds number 
flow past a thin wing at a moderate angle of attack, the viscous 
boundary layer and wake are essentially confined to a thin 
region. Outside the boundary layer and wake velocity gradients 
are small enough that the shear stresses acting on a fluid 
element can be neglected. For a region near the centre of 
rolled up region, however the velocity gradients are large and 
the viscous effects become dominant. In the present study only 
the initial roll up stage is analyzed and, consequently the 
fluid is considered to be inviscid. 
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The problem is fomula-ted relative to a wing fixed frame 
of reference xyz. The x~axis is the wing centreline and the 
xy plane is the wing plane of syrimetry V7ith z-axis in the 
spanwise direction (Fig. l). The dimension-less free stream 
velocity is expressed by 

e<» = cosa i -- sin a j (2.1) 

v/here i, j# k are the base unit vectors of the xyz frame of 
reference and a is the angle of attack. 

The continuity equation for an incompressible flow is 

V.V = 0 (2,2) 

and the irrotationality of flow demands 

T^x V = 0 (2-3) 

where V is the velocity vector. Equation (2.3) holds everywhere 
in the flow field except on the wake surface and on the lifting 
line where rotational regions are confined. If we define the 
velocity potential and perturbation velocity potential 4> 
such that 

V = Coo + V(|)= (2,4) 

Then Eq, '(2.3) is satisfied identically. Substituting Eq. (2,4) 
into (2,2), the Laplace equation for both ^ and 4* is obtained 

i# e . , 

v^<f= 0 * ^ (J) = 0 (2.5) 

The no-penetration on the wing surface S(r) = 0 relative 
to wing fixed frame of reference is given by 
(e« + V(|>). = 0 on S(r) = 0 


( 2 . 6 ) 
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the no*-»penet ration condition is applied to enforce the flow 
tangency on the vrakc surface 

(£oo + ^ $)• vf = 0 on f(r) = 0 (2.7) 

V7here f(r) = 0 is the equation of walce surface. 

The no-pressure jump condition on f (r) is obtained from 
incompressible Bernoulli’s equation 

cp(r) = -V({).[V(J)+ 2e«] ^2.8) 

where cp (r) is the pressure coefficient at any point (r) . 
Poming the pressure jump from Eq. (2.8) and etiuating the result 

to zero, we obtain 

Acp = cp, - op = -( V(l)^- ( V(^ + ^l |>2 + 2 £- ) = 0 

(2,8a) 


where subscripts 1 and 2 refers to upper and lower surfaces 
on the wing, respectively. 


Rearranging (2.8a) and setting 


•fl - <t>2 = 


(2.8b) 


Acp = - 2 (Vf .V)(A(j)) = -2 (A(l)) = 0 (2.8c) 

where is the velocity of wake element relative to xyz frame 
of reference 

Eq. (2.8c) represents theorem of Kelvin and Helmhaltz 
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5t ~ ” §t ~ * Ba ® f =0 (2*9) 

of conservation of circulation and outflow of vorticity, 
respectively. In Eq, (2.9) # n^ is a unit normal to the surface 
A bounded by a closed curve around which circulation is 
calculated. Eq. (2.9) simply states that the rate of change 
of circulation around a closed curve or the rate of change of 
outflow of vorticity through the surface bounded by this closed 
carve is zero (following the same fluid particles). 

For uniqueness of the solution one has to impose the Kutta 
condition along the edges of separation. Here Kutta condition 
is represented by 

Acp/ = 0 (2.10) 

TE 

finally the infinity condition requires that 

p* (J) — ► 0 away from f(r) and S(r) (2.11) 

The basic unknowns in the present problen are the vorticity 
distribution and free vortex sheet. They are determined by 
satisfying the Eqs, (2.1-2.11). Eq, (2.7) requires the flow to 
be tangent to wake surface while (2.8a) requires this tangential 
flow to be parallel to the vorticity direction. Therefore if 
position of f is adjusted so that flow direction is parallel 
to the vorticity direction on the surface f, the boundary 
conditions of Eq. (2.7) and (2.8a) are autcanatically satisfied. 



CHAPTER - 3 


METHOD CP SOLUTION 


3,1 General 

Uith the linear vortex lattice method (VLM) the so called 
horse-shoe elements are used to form the lattice as well as the 
lines modelling the %^ake adjoining the trailing edge. The directio 
of these lines is prescribed, not obtained as part of the solution 
end no attempt is made to model the tip or leading edge vortex 
system. Consequently the only unknovjns are the circulation of 
horse-shoe elements. These are found from the no-penetration 
boundary condition and consequertly the entire problem is 
linear. These methods yield accurate results for small angle 
of attack only. 

A steady non-linear hybrid vortex (NHV) method /SO/ is 
considered for large aspect ratio wings at large angles of 
attack. In this method vortex panels with linear vorticity 
distribution are used on the wing as well as in the near field 
of the wake. For the far field calculations, the distributed 
vorticity over each far field panel is lumped into equivalent 
concentrated vortex lines. In this way accuracy is satisfied 
in the near field, while computational efficiency is maintained 
in the far field. The coupling of a continuous vortex -sheet 
representation and a concentrated vortex line representation 
for solving the non-linear lifting surface probl^ is referred 
to as "non-linear hybrid vortex" (NHV) method. 
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The circulation distribution method /61/ is presented 
to calculate the wake roll-up for a steady flow past a high 
aspect ratio v/ing at lov7 angle of attack, in which the linear 
vorticity distribution on the wing is replaced by a quadratic 
circulation distribution over the lifting surface. The flow 
field considered here is incompressible and inviscid. The present 
wake model bears no small perturbation assumptions and all the 
equations engaged are solved analytically. Moreover the method 
is capable of calculating the wake geometry for a variety of 
planforms and shapes when the spanwise circulation distribution 
is known either fran nvunerical calculations or experiments (VliM 
is currently used) , As the nxamber of panels in the spanwise 
direction increases the v;ake sheet becomes more flexible and 
tends to roll-up more. To prevent chaotic motion at the tip, the 
method due to moore /50/ is used to model the tip core, 

3.2 Vorticity Distribution method 

A higher order panel method is used to calculate the three 
dimensional w^e roll-up behind a large aspect ratio wings in 
steady inccwnpressible flow. 

For a steady syrrsnetric flow the governing Eqs. are (2.1-2, 11) , 
The vorticity distribution over each bound vortex panel is 
obtained by satisfying the no-penetration condition over the 
control points, continuity condition at each common node of 
bound vortex panel, Kutta condition at the trailing edge and 
symmetry condition along the root chord assuming wake to be 
planar (Fig. 2) • 
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Equating the outflov/ of vorticity from bound vortex 
panels to the free vortex panels along the separation edges / 
induced velocity at each corrmon nodal point over the wake is 
calculated in the near region which is taken two spans in length. 
In the far field the distributed vorticity is reduced to 
concentrated vortex lines and simple Biot-Savart lax-r is used to 
calculate induced velocity due to vortex filaments, see 
Appendix A, Wake geomet3ry is then evaluated in an iterative 
fashion by satisfying conservation of circulation, flow tangency 
and zero pressure jump for all v;ake panel surfaces in the near 
region. 

Planar quadrilateral vortex panels are used to represent 
the bound vortex sheet. On each panel a local linear distribution 
of vorticity in term of five undetermined coefficients (a^ “ 
is specified in local axes systen, in vAiich S and I axes are 
located in the panel plane such that the 5-axis coincides with 
the side 1-4 of the quadrilateral panel. The ti axis is perpendi- 
cular to the plane such that I and tj form a right handed 
local coordinate system. See Pig* 3a, 

Wg = a^ + Sa 5 + 33 I (3.1) 

W| = “'^2^"’ ^4 ^5^ (3.2) 

This type of yorticity distribution satisfies solenoidal 
property of vorticity 


y.w = 0 


(3.3) 
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Planar triangular panels are used to model free vortex 
sheet where highly nonplanar and twisted surfaces are encountered. 
The local vorticity distribution is still given by Eqs, (3.1) 
and (3.2). The comers of this panel serves as nodal point. 

The ?-axis now coincides with side 1-3 of the panel, the 
I -axis is in the panel plane and n axis is perpendicular to the 
plane such that I, and t) axes again form a right handed 
local coordinate system. See Pig. 3b, 

The velocity induced by a distribution of vorticity on 
the surface a at any point P(x,y, z) in the local panel coordinates 
can be calculated by means of Biot-Savart law to yield 

V (x,y, 2 ) = - 7 — — ZL_ aa (3.4) 

(r-s)^ 

where r = x e|+ ye^ + ze ^ is the position vector of a point P, anc 
s =w I© ^+w^e ^position vector of a point on q. With the vorticity 
distribution 


w = wtefc + wye 


(3.5) 


the induced velocity in the near field at any point is given by 

^ywge|+ ((z-l)w| - (x-l)wg ) e^-yw | e ^ ^ ^ 

((X- 1)2 + y2 + (2-5)2)3/2 

Each rectangular panel over the wake surface is divided 


V (x,y,z) 


into two triangular panels. Assuming the same vorticity 
distribution over both the panels (Pig, 1), Eq* (3,6) can be 
expressed as 
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V (x,y, 2 ) 


1 

4 % 


2 

Z 

i=l 


ywg ®| -(x-l)w 

/j ^ 

((x-l )2 4 . y 2 ^j2j3/2 


(3.7) 

Substituting Eqs, (3,1) and (3.2) into (3,7) the induced 
velocities are calculated in local axes system in terms of the 
coefficients of a‘s# see Appendix 3. At each control point of 
bound vortex panel no-penetration condition Eq. (2.6) is enforced 
(Appendix B and C) in terms of five unknowns per panel. Control 
point is defined as the average point. Additional equations are 
obtained by satisfying conditions of vorticity continuity, Kutta 
and symmetry conditions, see Appendix C, 

The number of equations from these four conditions is given 
by 5NB + 2NZX(15X-1> + 2NXX{NZ-1) in term of 5NB unknowns, where 
NB is number of bound vortex panels. The vorticity distribution 
coefficients over bound vortex panels are obtained by solving this 
overdetermined set of equations using method of Least-Squares. 
Vorticity distribution coefficients over wake surface are derived 
in term of vorticity distribution coefficients of bound vortex 
panels, see Appendix D. 

3.3 Circulation Distribution Method 

To calculate the three dimensional wake geometry in a 
steady incompressible flow behind a symmetric wing of large 
aspect ratio for a given spanwise circulation, method uses 
triangular vortex panels with a linear vorticity distribution 
for a near wake modelling. To model this problem the wake is 
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first initialized as a flat sheet of zero thickness/ extended 
from the wing TE v;hich is divided into three regions-ad joining# 
near and for region. 

The adjoining region is described as a small region 
(l/lO of span) of planar vortex sheet where the wake surface 
lies in the direction of x-axis and stays fixed. The region 
is introduced to model the satisfaction of Kutta condition, 
that the flow will leave the trailing edge smoothly. 

In the near region approximately two spans in length, the 
wake surface is split into a netv?ork of flat, triangular 
elements (panels) with linearly distributed vorticity in local 
spanwise direction and constant vorticity along the strearawise 
direction (local flow direction) for each panel. 

In the far region, straight semi-infinite vortex filaments 
of equivalent strength are shed from the trailing edge of vortex 
panel sheet to infinity. Such representation for the far wake 
is introduced for economic feasibility while accuracy is 
maintained in the near region as well. 

The wake geometry is then evaluated in an Iterative fashion 
by satisfying conservation of circulation, the flow tangency 
condition and zero pressure jtunp for all wake panel surfaces in 
the near region. The iteration continues until the wake shape 
converges or equivalently the wake sheet truly becomes a stream 
surface in the near region. Straight semi-infinite vortex 
filaments in the far wake are assumed to be in the uniform 


stream direction 
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The governing Eqs. are again (2.1-2.11). The spanwise 
circulation distribution P(z) over the lifting surface is 
obtained by VLM approach. Once tlie spanv/ise circulation Piz) 
is obtained the vorticity shed by the V7ing or the strength of 
the trailing vortex sheet can be expressed as 


w= - 


dz 


(3.8) 


The direction of shed vorticity is taken perpendicular to 
lifting line. 


For a piecewise linear distribution of vorticity in the 
spanwise direction# along with the fact that die change of 


circulation { 


must be shed behind and into the wake. 


Equation (3.8) can be written with the use of trapezoidal arule 


+ ( 


n 






a ) 


( 3 . 9 ) 


The value of vorticity at each nodal point can be 
obtained by starting fron the symmetry plane and marching in 
the spanwise direction with the fact that 

w^ = 0 at z = 0 { i=l ) (3,10) 

and referring to Pig, 4, Eq. (3,9) reflects that the vorticity 
distribution is adjusted such that the shaded area under the 
vorticity curve is always equal to change of circulation 
provided that the vorticity is continuous at the ccwroon node 
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between two adjecent panels. Once the vorticity at each of 
the common node is knovm, it can be expressed as 

w = (a + B 5 ) e^ 


(3.11) 


The velocity induced by a distribution of vorticity 
on the surface a at the point P y^^# in the local panel 

coordinates can be calculated by using 2q. (3.4), Fig. 5. 

) (y^ e y - (z^- ^ ^ ^Ti ^ ^ 




0 ((x„-^)^ * 


(3.12) 


with substitution 

(x.-n = (y^ + tane 


(3.13) 


K tan © 


Equation (3.12) on integration gives 


5 (^) 

TE 


(A+Bl ) (Y™ (z -5)e ) 


1 r r - ^ m — i 

4ic q 2// t\2 2 / n2\% 

0 ^2 ((x^- 1)^ + y^ + (2^- a ) ^ (I) 

LE 


3,4 Free Vortex wake 

Ihe induced velocity at a point on the wake surface is 
calculated by summing the contribution from three major sources 

They are 

1) Bound vortex segment (Appendix A) 

2) Triangular panels in the near region (Sec, 3.2 & 3.3) 

3) straight semi-infinite vortex filaments in the far region 

(Aj^endix A) 
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The induced velocities calculated due to near wake 
region from Eqs. (3.7) and (3*14) have to be transformed to 


global system of reference as follov/s 




(3.15) 


where d with suffix is the cosine of angle between its 
two subscripts. 

To satisfy the boundary conditions over free vortex panels 

an iterative technique is followed. During a typical iterative 

cycle vorticity distribution is calculated. This is followed 

by aligning the surface, after perturbation velocities are 

calculated at the nodal point 1 and control point of panel for 

vorticity distribution and circulation distribution methods 

respectively. For calculating the induced velocity at the 

—6 

point results are evaluated at a distance of 10 

Each of the triangular panel is considered to be a stream 
surface ccanposed of the streamlines that are parallel to one 
another with respect to side 1-3^ Once the induced velocity at 
a nodal point or control point is determined the streamline 
iflhich contains that point is moved in the same direction as that 
of local velocity. 
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To demonstrate the approach, the unit vector of streamline 


at a control point, e is given by 

*"'**C]^ 


= = H. IS . ^ ^ 

-=p Vj, y 


e. 


(3.16) 


VJhere Uj^, v^, are the ccmponents of induced velocity in 

x,y & 2 direction respectively and V = (u^ + v„ + w„ ) ^ 

K K, K K 

Side 1-3 is then moved in the direction parallel to unit vector. 

Up 

(L+1, J) = x^ (L, J) = x^ (L J) + 113 ^ 3 ! ~ 

V 

y^ (L+1, J) = y 3 (L, J) = y^ (L, J) + II 13 I — 

w 

2 ^ (L+1, J) = Z 3 (L, J) = 2 ^ (L, J) + ll^^j ~ 

^R 

VJhere (L,J) and (L+1, J) refers to upstream and downstream 
nodes of side 1-3, 1^^^ is the length of side 1-3 (Pig. 5) • 

The best approximation for the angle at which flow leaves 
the trailing edge is half the angle of attack. The wake is 
over-released by taking the first row of panels in the direction 
of chord for first iteration. The computation is initiated by 
specifying the rest of the wake to be planar. Align the second 
row of panel (Pig. 6 ) with local velocities. The downstream wake 
shape is adjusted so that the nodal points will have the same y 
and z coordinate as that of second row. The vorticity distribution 
is calculated again due to new wake geometry. Align the next row 
with local velocities and process continues until last row of 
panels is moved, ihe semi- infinite filaments of far wake are left 
parallel to free stream velocity. 



starting from the second iteration the process is 
altered due to the fact that the geometry of the upstream wake 
will be affected. Calculate the resultant velocities at all of 
the nodal points and move the dovnstream nodal points simulta- 
neously in the direction of local velocity. Repeat the same 
process until v/ake shape is corverged. 


3,5 Tip Core modelling 

Refinement of the wake roll up /50/ is done by dividing 
the spiral into inner and outer portions. Vortex sheet is 
represented by NF number of points in the spnnwise direction. 
We assume that a turn of spiral is represented by M point 
vortices. This N has to be determined by trial and error, 
actually N/NP should be 0,1 to 0,15, 


Let e 


360 

N 


(3,17) 


We first calculate induced velocities and new coorotinates 
of y and z, and then angle 6betv?een the two vectors of adjacent 
panels is determined. If Q exceeds 0^ then two panels are 
combined to form a new tip vortex placed at the centroid of 
those two vortices with circulation that equals the combined 
circulation of original pair. 

We start our calculations from the wing root. Initially 

c which is called outer part of the spiral 

and takes about 20% of the vortices. When 6>0 then all the 

c 

vortices are combined to form a tip vortex or inner portion of 
the spiral, which also constmies almost the same amount of vorticei 
©= 0^, can be Included in either inner part of the spiral or 
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outer part. It should be noted that once the condition 0>© 
is achieved then for all the points on the right side, the same 
condition has to be satisfied. Vortex sheet will came nearer 
and nearer to the tip vortex as go on increasing ::P, 


Referring to Pig, 11, let 


'i = h-i - = ‘^i-i - 


= Y-i j + ^ 

1 X 


where is a two dimensional vector which represent the 

coordinates (y(i), z{i) 

and 


^i+1 


= X 


X 


i+1 


^i+l J 


-f k 


the angle between the two vectors and is given by 

e = cos-1 yj yi-n ^1 

(y| + + 4-n^^ 

whose direction cosines are ^y. , z . > and <y. z.,.>. 

X X x-rl X+1*^ • 



CHAPTER - 4 


DESCRIPTION OF COKPUTBR PxROGRAM 

A computer program is developed to implement the method 
of solution for both the methods. Vorticity distribution 
method has been used to calculate the normal force coefficient 
over flat plate at an angle of attack assiaming wake to be 
planar. The wake rollup could not be solved because of the 
limitation on the computing pov/er available. The computer 
pix>gram consists of several function subprograms and subroutines. 

The subroutine CORT calculates the nodal point of panels 
over the flat plate wing with given aspect ratio and number of 
panels in both the direction as input. 

Induced velocities are calculated to apply the no-penetratioi] 
condition over the average point of each paiSel. This has been 
achieved by using f traction subprograms SIMPSN and CINT'S, Next 
letter of CINT i,e, 5th denotes the j corresponding to a^ such 
as CINTIR will give the y— component of induced velocity per 
unit a^ for right wing and CINTIL for left wing. SIMPSN integrate^ 
all these functions by Simpson rule. 

Local coordinate axis system is needed for the calculation 
of perturbation velocities, otherwise it is difficult to 
integrate those expressions of induced velocity in global syston 
of reference. This way lower and upper limit for integrals of 
induced velocity will be same for all the panels, provided the 
size of panel over vhich integral is taken is same. By the 
shifting of the origin at one of the nodal point (comer point) 
coordinates of panels in their own axis system are calculated 
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and then used for tha calculation of induced velocities. 

These induced velocities are trans formed bach to the vring 
fijied frame of reference by rotating the axis system. The 
induced velocities calculated do not require any transformation 
as the v;ing taken is a flat plate, 

CEITIR to CIIITSR are used for the calculation of induced 
velocity in y-direction and CIIIT6R to CIImTSR are used to 
calculate the induced velocities in x and z-directions. 

Function subprograms ZO and ZU are the equations of line which 
have been used to convert quadrilateral panels to triangular 
panel with local origin at the left and right mode respectively. 

Continuity of vorticity in x and z directions, Kutta 
condition along the trailing edge and symmetry condition for 
a symnetric flow along the root chord are applied, details of 
which are given in Appendix - C, All conditions including 
no-penetration condition have been stored in the matrix A<in,n), 
where m is the total number of equations in term of n( 5 times 
the number of bound vortex panels) variables B is a column 
vector. This overdetermined set of equations (since number 
of equations are more than the ntamber of variables) can be 
solved using the method of Least-squares by reducing them to 
their normal equations. This schene was not acceptable to the 
DEC-1090 system for some cases of increased number of panels 
as it required much more memory space. Therefore standard 
library subroutine I MSI* LLSQAR was selected to solve set of 



ecjuations by the method of Leas t-Srraa res . Nov? the vorticity 
distribution coefficients so obrained are stored in matrix 3 
v;hich over\-?rites the colvimn vector of right hand side. This 
problem is solved assuming wake to be planar. 

Ccxnputer program for the circulation distribution method 
is divided into three major parts and each part consists of 
several subroutines and function subprograms. 

The first part deals with the v/ing geometry, wing 
panelling (bound vortex panels), wake panelling (free vortex 
panels) and calculation of circulation distribution over the 
lifting surface. The various input to this part are aspect 
ratio, angle of attack, number of panels in the spanwise 
direction of flow over lifting surface and free surface. 

Subroutine CORD calculates the coordinates of each panel 
and control point over the lifting surface in order to calculate 
circulation distribution, while coordinates of triangular panel 
over the force vortex surface are determined by the subroutine 
GEMTRY assuming wake to be planar. Subroutine GMQUAD serves the 
same purpose but the length of panels in the streamwise direction 
is varied according to quadratic law. Panels in the direction of 
span are distributed according to cosine spacing for both subrou- 
tines GEMTRY and GMQUAD. 

No-penetration condition of tangential flow is applied at 
the control point of each panel over the flat plate. Coefficient 
matrix (w) obtained frc»n this condition is solved by subroutine 
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Flirv which has been tested against the standard NAG subroutine# 
FOIAAF. Both of these (PClAr'vF and FIlJV) calculate the inverse 
(UNIT) of matrix (w) and this UNIT is multiplied to the right 
hand side of simultaneous equations to obtain the unknovms 
of circulation. 

Second part deals V7ith the v/ake. The wake geometry is 
generated by dividing the straight xvake into NVJ rows and RZ 
columns (Semispan) to form (NW) (RZ) rectangular elements {Pig-6) , 
Each rectangular element is then divided into tv70 triangular 
panels. Once the spanwise circulation is obtained, vorticity 
distribution over each triangular panel is calculated. Function 
subprogram SIMPSN integrates CINTIL and CINTIR to give y-ccanponent 
of induced velocity due to left and right wing and corresponding 
2 -ccanponents are obtained frcra CINT2L and CINT2R. To calculate 
the velocities in the global axis system, induced velocities 
calculated are multiplied by the transformation matrix. Induced 
velocity at a point over the wake surface is the sum of contri- 
bution f rxxn bound vortex panels , far wake region and near wake 
with triangular vortex panels. Contribution due to bound and 
far wake is calculated by VELY, VELZ, VLINFY and VLINFZ, Function 
subprogram VELY and VELZ will give induced velocities in y and 
2 -directions due to finite vortex element with given end points 
while VLINFY and VLINFZ calculate velocities due to semi-infinite 
vortex filament in y and 2 -directions respectively. Aligning 
the wake surface in the direction of local flow velocity is 
done itself in the main program. 
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In the third part of program nxmiber of vortices per 
turn of rollup is tested. If number of vortices per turn 
falls from some specified level then those vortices are 
combined to form a tip vortex and coordinates of other 
vortices which do not satisfy this condition are written as 
it is. 
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VARI/\BLE 

A 


ALPHA 

B 

D 


EXL 

P 


G 


IDGT 

lER 

IT 

ITER 

N 

NEQ 

MJNK 

m 

NX 

NP 


Vorticity Distribution Method 

DESCRIPTION 

mxn matrix, with m number of equations in 

term of n unknoms 

angle of attack in degrees 

col-umn vector (See IMSL : LLSQAR) 

is a 3x1 matrix for the three components 

of unit normal at the receiver panel 

value of J 

is 3x3 matrix obtained from the dot product 
of global and local axis system 
is N X NUNK matrix for the no-penetration 
condition over the average point of bound 
vortex panels 

nimber of significant digits for IMSL 
error parameter 

do loop index, performing iterations 
maximum number of iterations allowed 
total bound vortex panels 
total number of equations 
number of unknowns 

panels in the direction of wake over the 
wake surface 

number of panels in the direction of flow 
over flat plate 

number of nodal points per unit s«ni-span 
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VARIABLE 

NZ 

512 

513 
VRl 


VLl 


XP 

XAyl 

XAy2 

XAy3 

XX 

yyL 


yyR 

ZZl* 

zzaR 

zeta 

ZP 


DESCRIPTION 

semi-spanwise niomber of panels ; 

distance between the points 1 and 2 
distance between the points 1 and 3 
a coltamn vector v;ith arrays equal to nximber 
of bound vortex panels and calculates induced 
velocity in y-direction for the right hand 
side of wing per unit a^ 

- similar to VRl but for left wing - 

induced velocity per unit a^ in local coordinate 

system 

width of panel in x-direction 

XP for first row over the wake (adjoining region 
XP for 2^*^ row 
XP for 3^"^ row 

x-coordinate in local axis system for left 
as well as ri^t wing 

y-coordinate for left wing# vhen observer 
is facing the flow 

y-coordinate for ritjht wing in local axis system; 
z-coordinate for left wing 
z-coordinate for right wing 
value of 

width of panel in z-direction 
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Circulation Distribution Method 


VARIABLE 

A 

AA 


ALPliA 

AR 

AREA 

B 

BB 

C 


CR 

CT 

E2 

E3 

E4 

P*s 


(M 

GH 


DESCRIPTION 

lov/er limit of integeration 

constant term of the linearly varying 

vorticity expression 

angle of attach in degrees 

aspect ratio 

area of half flat plate 

upper limit of in teg ration 

coefficient of linearly varying vorticity term 
It is a column vector. Right hand side of 
simultaneous equations obtained from no— penetr< 
ion condition 
root chord 
tip chord 

widtb of panel in spanwise direction 

margin left around the singularity on one side 

margin left on the right side of singularity 

F is a 3x3 matrix obtained from the dot 

product of local and global axis and elCTients 

are like P(l,l) as Pll 

circulation obtained by joining the panels 

in the direction of flow 

circulation over each panel 

row number for free vortex panels 


II 
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VARIABLE 

IT 

ITER 

NB 

NF 

NP 

NW 

NX 


Ny 

S13 

THETA 


THETAc 


TU 


TV 


TW 

UNIT 


VB 


DISCRIPTION 
iteration n\mber 

maxiimim number of iterations allowed 

total number of bound vortex panels 

free vortex panels 

total nodal points per serai-span 

number of panels in the direction of flow 

number of panels in the flow direction over 

lifting surface 

spanwise number of panels 

distance between^ the points 1 and 3 

angle between the line joining both sides of 

nodal point in order to obtain the smooth 

roll-up 

critical angle, if TOETA exceeds THETAc 
vortices are combined 

total velocity at the control point of a 
triangular panel in x-direction 
total velocity in y-direction at the control 
point of a triangular panel 
total velocity in z-direction 
inverse of matrix w obtained from no- 
penetration condition over each panel of 
lifting surface 

y-ccanponent of induced velocity due to 
bound vortex panels 
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VARIABLE 


DISCRIPTION 


VF 


vw 

VLy 


VLz 

VRy 

VRz 

W 

WF 

WW 

XP 

XX 

yyL 

yyR 

EEL 

ZZR 


y-component of induced velocity due to far wake 
( semi-inf inite vortices ) 
y-c<xnponent due to near wake region 
induced velocity in y-direction due to left 
wing in local axis system at one control point 
due to single triangular panel and obtained by 
integ- rating function CINIIL 

- z-di recti on due to left wing - 

- y-direction due to right wing - 

- z-direction due to right wing - 
z-component of induced velocity due to bound 
vortex panels 

z-component of induced velocity at a triangular 
panel due to all the far wake vortices 

- due to near wake - 

width of each panel in x-direction 
x-coordinate in local axis syst^i 
y-coordinate in local axis system for left wing 
y-coordinate in local axis system for right wing 
z-coordinate for left wing 
z-coordinate for right wing 



CHAPTER -» 5 


RESULTS AND DISCUSSION 

The methods of vorticity distribution and circulation distributior 
are implemented using a computer program. The wake shape from 
vorticity distribution method could not be obtained as the 
necessary computing power was not available. Spanwise variation 
of pressure coefficients over flat plate rectangular wing of 
aspect ratio 1.0 at an angle of attack 15° are obtained and 
presented in Pig. 12. The results are compared with those of 
Kandil et al /52/ who have considered the wake roll up. 

The difference between the two curves at a section 
perpendicular to the free stream direction, as shown in Pig. 12 
is due to the planar wake for our analysis. 

In the program for circulation distribution method the 
initial wake geometry is autcroa tic ally generated by dividing 
the straight wake into NW rows and Nz columns to form (NwXNz) 
rectangular elements. Each rectangular element is then divided 
into two triangular panels. The panels in the near wake are 
distributed using cosine spacing in the spanwise direction to 
emphasize the details near the wing tip. For all the numerical 
runs presented, the results of the wake shape are shown only 
for the right half of the wing, since the flow considered is 
symmetric , 

The convergence properties of the circulation distribution 
method are studied for a rectangular flat plate wingcf AR=8.0 at 
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5-deg* angle of attack and panels in the spanwise direction 
are varied from 6 to 18 as shovm in Pig. 13. It is seen that in 
the third iteration the nodes of the wake shape at 1.69 span 
lengths behind the trailing edge virtually fall on top of one 
another and form a defined contour. The wake geometry shows 
divergence after three iterations in each case. Such divergence 
is due to truncation errors. This has been checked by reversing 
the direction of vortices after three iterations and original 
wake sheet was not produced. 

It is also learned that as the number of spanwise panels 
increase bringing the outboard panel closer to the tip the wake 
sheet becomes more flexible and tends to roll-up more. This is 
because the present method is limited to Calculations with less 
than one turn of rollup. 

Simpson rule is used to evaluate the integral, Eq. (3.14), 
for the expression of induced velocity. Singularity contributes 
to zero induced velocity, so we can leave an ecpaal margin on 
both the sides of singularity to calculate the contribution of 
the panel towards induced velocity having singular effect. The 
results of Fig. 13 are computed by leaving an interval of 1/30 
of panel width (spanwise) . The results are ccmnputed by selecting 
two different intervals 1/30^ and 1/50 of panel width and 
are plotted in Pig. 14. It is seen that the wake shape is not 
changed by selecting different intervals, however the panel end 
points near the wake tip show different magnitude. 
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For the above and the follo\:ing calculations the length 
of the near wake taken is txra spans. Pig. 15 illustrates the 
pattern of rollcd-up wake sheet in the streamwise direction for 
a rectangular wing of AR = 8.0 at 5-deg. angle of attack, for 

five different dovmstream stations, x = 0.29, 0,59, 0.89, 1.19 

w 

and 1,49 span lengths measured from trailing edge; panel 
arrangement of Nx = 1, Nz = 12, Nw = 7 is used. The streching 
and spiral motions are seen on the v;ake panels over the tip region 
during the roll-up process xvhile the v/ake sheet remains flat at 
the middle of the wing, but y-coordinate increase in the streamwisi 
direction. 

The effect of changing Nx, the nvimber of panels over the 
lifting surface (used to calculate the circulation distribution) 
in the streamwise direction, over V7ake geometry is negligible 
as shown in Fig, 16, 

The results of wake shape for flat plate rectangular wing 
of AR = 8,0 at two different angles of attack, iK = 5 and 10 deg, 
are shown in Fig. 17. A comparison of the wake roll up with the 
results of Sucio and Morino /51/ at 5 chord lengths downstream 
of the trailing edge, demonstrates the validity in predicting 
the wake roll-up. The agreement is satisfactory in the sense that 
Sucio and Morino used constant doublet panels to model the wake; 
this is equivalent to using vortex filaments, Sucio and Morino ’s 
result for the location of the vortex centres appear to agree 
well with experimental observations oc Chigier and Corsiglia /40/ 



For spanwise panelling Kz = 6-18 the near wake exhibit no 
more than one turn of roll-up at any streamwise location. Pollotij-i 
Moore /50/ we tried to model the core by increasing the number 
of panels in the spam-.’ise direction. Increasing nxmber of panels 
in the spanwise direction beyond 20 results in overlapping of 
panels or in other words jumping of panels over one another 
results in the tip region and no definite contour was obtained. 



48 


LIST OP RSPER.'.NCES 


1, Vom T, Karman (1911) Hydrodimamics by H, Lamb (1945) Dover 
Publication P-225, 

2, Prandtl (1918) Principles of Ideal-fluid aerodynamics by 
Karamcheti, K, (1966) John Wiley & Sons, Inc,, N.Y. 

3, Roshenhead, L. (1931) The formation of vortices frcan a 
surface of discontinuity, Proc. R,Soc. bond. A-134,pp, 170-192, 

4, Westwater, P,L, (1935) Rolling up of the surface of 
discontinuity behind an aerofoil of finite span, ARC R&M 
no, 1692, 

5, Bollay, W, (1937) A theory for rectangular wings of small 
aspect ratio. Journal of Aeronautical Sciences, vol, 4, 

294-296, 

6, Multhopp, H, (1950) Methods for calculating the lift 
distribution of wings (subsonic lifting-surface theory) 

ARC R&M NO 2884, 

7, Lawrence, H, (1951) The lift distribution on low aspect 
ratio wings at subsonic speeds. Journal of the Aeronautical 
Sciences, Vol, 18, 683-695, 

8, Brown, C,E, and Michael, W,H, (1954) Effect of leading edge 
separation on the lift of a delta wing. Journal of Aeronauticaj 
Sciences, Vol. 21, 690-694. 

9, Mangier, K,W. and Smith, J.H.B. (1957) Calculation of the 
flow past slender delta wings with leading edge separation, 

RAE Report No 2593, 

10, Gersten, K. (1961) Calculation of non-linear aerodynamic 
stability derivatives of aeroplanes, AGARD Report No 342, 

11, Bir)choff, G, (1962) Helmholtz and Taylor instability 
Proceedings of the Symposium on Applied Mathematics, Vol. 13, 
55-76, American Hathafnatical Society. 

12, Rubbert, P.E. (1962) Theoretical characteristics of arbitrary 
wings by a nonplanar Vortex Lattice Method. The Boing Company 
Report D6-9244, 

. Gamer, H.C, & Lehrian (1963) Nonlinear theory of steady 
forces on wings with leading-edge flow separation, ARC R&M 
No 3375. 


13 



49 


14. Van Dyke, M, (1964) Liftin'^ line theory as a singular 
perturbation problem. Arch. Kech. Stos., Vol. 16, No. 3. 

15. Takami, H. (1964) A numerical experiment v/ith discrete 
vortex approximation with reference to the rolling up of a 
vortex sheet. Department of Aeronautics and Astronomy, 

Stanford University, report SUDAER 202, 

16. Thurber, J.K, (1965) An asymptotic method for determining 
the lift distribution of a s%rept~back v/ing of finite span, 
comm. Pure Appl, Math. 

17. Ashley, H, and Landahl, M, (1965) Aerodp:-.=?r.ics of wings and 
bodies. Addison-Wesley Publishing Conpany I.C,, Mass. 

18. Kuchemann, D. (1965) Report; on the I.U.T.A.M. symposium on 
concentrated vortex motions in fluids. Journal of Fluid Mech 
Vol. 21, 1-20, 

19. Ermolenko, S.D, (1966) Nonlinear theory of small aspect 
ratio wings. Soviet Aeronautics (in English) Vol. 9, 5-11. 

20. Hedman, S.G. (1966) Vortex Lattice Method for calculation of 
quasi steady state loadings on thin elastic wings. The 
Aeronautical Research Institute of Sweden, Report 105, 
Stockholm, 

21. Belotserkovskii, S.M. (1967), Theory of thin wings in subsonic 
flow. Translated from Russian, Plenxam Press, N.Y., 67, 

22. Kalman, T.P,, Rodden, W.P. and Giesing, J.P. (1967) Aerodynamid 

influence coefficients by the doublet lattice method for | 

interfering non planar lifting surfaces oscillating in 
subsonic flow (Part I) , McDonnell Doughlas Corporation, 

Report No. DAC-67977. 

23. Sack, A.H, , Lundberg, R.E., and Hanson, C.W. (1967) A ; 

theoretical investigation of the aerodynamics of slender wing : 
body ccHnbinations exhibiting leading-edge separation NASA : 

CR-719. 

24. Landahl, M.T,, and Stark, V.J.E. (1968) Numerical lifting- 
surface theory problems and progress, AIAA Journal, Vol, 6, 
2049-2060. 

25. Belotserkovskill, S.M. (1968) Calculation of the flow around 
wings of arbitrary planforms over a wide range of angles of 
attack. NASA TT P-12, 291, May 1969. 

26. Smith, J.H.B, (1968) Improved calculations of the leading 
edge separation from slender, thin, delta wings, Proc, 

R.Soc, Lond,, A 306, 67-90. 

27. Albano, E,, and Rodden, W.P, (1969), A doublet-lattice method 
for calculating lift distributions on oscillating surfaces in 
subsonic flows. AIAA Journal, Vol. 7, 279-285. 



50 


28. 


29. 


30. 


31 


« 


32. 

33. 

34. 

35. 

36. 


37. 

38. 

39. 

40. 


uames# R.M. (1969) On the remarkable accuracy of the vortex 
latt5.ce discretisation in tliin v;ing tiieory. KcDonnell Doughl 
7*ircraft Co., Report i)AC-67211, 


P n. 


Cunnighara, ‘-i.K, (1971) kn efficient subsonic callocaticn metho 

for solvinc lifting-surface oroblemG. Journal of 7-i.ircraft, 

Vol. 8, 168-176. 


Polhamus, 3.C, (1971) Charts for predicting tns subsonic 

vortex-lJft ch3iracteristics of Arrov;, Delta & Disiriond u'ings 
I-TASA TI'JD - 6243. 


Folhamus, E.C, (1971) Prediction of vortex-lift characteristic 
by a loading edge suction analogy. Journal of Aircraft, 

Vol, 8, 193-199. 

iioore, D.W. (1971) The discrete vortex approximation of a 
vortex sheet. Airforce office of Scientific Research Report 
1084—69, California Institute of Technology. 

Butter, D.J. and Hancock, G, J, (1971) A numerical metliod for 
calculating the trailing vorte:: system behind a swept wing 
at low speeds. The Aeronautical Journal, Vol, 75, 564-56S, 

Ashley, H. and Rodden, V7,P, (1972) IJing-body aerodynamic 
interaction. Annual review of Fluid Mechanics, Vol, 4, 431-472, 

Finkleman, D, (1972) Nonlinear vortex interaction on taring- 
canard configurations. Journal of Aircraft, Vol. 9, 399-406, 

Giesing, J.P., Kalman, T.P, and Rodden, VJ.P. (1972) subsonic 
steady and oscillatory aerodynamics for multiple interferring 
wings and bodies. Journal of Aircraft, Vol. 9, 693-702. 

Kentey, K.P, (1972) A correction to lilting line theory as 
a singular perturbation problem, AIA/i Journal, Vol. 10, 1683— 16S 

Rodden, W.P,, Giesing, J.P, and Kalman, T.P. (1972) Refinement 
of the nonplanar aspects of the subsonic doublet-lattice liftif 
surface method. Journal of Aircraft Vol. 9, 69-73. 

Hess, J.L. (1972) calculation of potential flow about 
arbitrary three dimensional lifting bodies. McDonnel Douglas 
Corporation, Report No MDC J5679— 01, 

Chigier, N.A. and Corsiglia, V.R, (1972) wind-Tunnel studies 
of wing wake turbulence. Journal of Aircraft Vol, 9,820-825, 


41 


Hua, H.M. (1973) A finite elenent method for calculating 
aerodynamic coefficients of a subsonic airplane. Journal of 
Aircraft, Vol, 10, 422-426. 



51 


42. 


43. 

44. 

45. 


46, 


47. 


48. 


49. 


50. 

51, 


52. 

53. 

54. 


1 !. T.. KA^iPUR 


ACC. Ko. 


■Jijiir 


/o 


Rom, J,, and Zorea, C, (1973) The cslcula-tion of the lift 
distribution and the near ■^rahe behind hiqh and low aspect 
ratio wings in subsonic flo;,’, Technion- Israel Institute of 
Technology, Haifa, 2AS Report 168. 


Kuo, C.C,, and Morino# L. (1973) steady subsonic flcv/ 
around finite thickness vrincs, Boston University, Department 
of Aerospace Engineering, Boston, Mass., TR-73-62. 

Kossow, V.J, (1973) On the inviscid roll-up stiructure of 
lift generated vortices. Joum-ial of Aircraft, Vol. 10, 647-65 


Bradley, R.G. Smith, C.W. and Ehateley, I.C. (1973) Vortex- 
lift prediction for complex wing planfoxms. Journal of 
Axircraft, Vol, 10, 379-381, 


Hough, G.R, (1973) Remarks on vorte2<-lattice method. Journal 
of Aircraft, Vol, 10, 314-317. 


Morino, L. and Kuo, C.C, (1974) subsonic potential aerody- 
namics for complex configurations : A general theory, AIA/v 
Journal, Vol. 12, 191-197. 

Mook, D.T, and Maddox, S.A. (1974) Extension of a vortex 
lattice method to include the effects of leading-edge 
separation. Journal of Aircraft, Vol, 11, 127-128, 

Rodeen, W.P,, Giesing, J.P., Kalman, T.P., and Rowan, J.C, 
(1974) Comment on : A finite element method for calculating 
aerodynamic coefficients of a subsonic airplane. Journal of 
A.ircraft, Vol, 11, 366-367. 

Moore, D.W, (1974) A numerical study of the roll up of a 
finite vortex sheet, J, Fluid Mech,, Vol, 63, 225-235, 

Sucio, E.O, and Morino, L. (1976), A nonlinear Finite 
element analysis of wings in steady incompressible flows 
with wake Roll up. AIAA paper 76-64, 1976. 

Kandil, O.A,, Mook, D.T., and Hayfeh, A.H, (1976) Nonlinear 
Prediction of the Aerodynamic loads on lifting surfaces. 
Journal of Aircraft, Vol. 13, 22-28, 

Pink, P.T. and Soh, W.K. (1978) A new approach to roll up 
calculations of vortex sheets, Proc, Roy. Soc, of London, 
Ser A 362, 195-209. 

Johnson, P.T., Tinoco, F.N., Lu, P., and Epton, MA. (1980) 
Three dimensional flow over wings with leading edge vortex 
separation, AIAA Journal, Vol, 18, 367-380, 



52 


55, Lavin, D, and. Katz, J, (19£1) vortex lattice method for 
the calculation of the ncnsteady separated flo\; over delta 
wings. Journal of Aircraft, Vol. 15, 1032-1037. 

56, Ktxman, S.M, and Stremel, ?,K. (1987) A vorteic wahe caoturing 
method for potential flow calcularions. AI.Al Paper 82-0947,193 

57, Svimma, J.M, (1982) Advanced Rotor analysis method for the 
aerodynamics of vortex/blaoe interactions in Hover. Eighth 
European Porum, paper - 23, Prance. 

58, Hoeijmaker, H.W.M, (1983) Computational vortex flow A.erodynami 
AGARD CP 342. 

59, Hoeijmaker, H.W.M. and Vaatstra, W*. (1983) h higher order pane 
method applied to vortex sheet roll up, AIA.-. Journal, Vol. 21, 
516-523. 

60, Kandil, O.A, , chu, L,, and Tureaud, T,, (1984) A nonlinear 
hybrid voirtex method for xjings at large angles of attack, 

AIAA Journal Vol. 22, 329-336, - 

61, Yeh, D.T, , and Plotkin, A, (1986) Wake roll up behind a 
large aspect ratio wing, AliV* Journal, Vol. 24, 1417-1423, 



53 


APPEKDK - A 


Calculation of induced velocity due to vortex filaments 

The far wake is modelled by straight semi-infinite 
vortex filaments that are in the uniform stream direction. 
Circulation p for this is derived in terms of the unknovm 
parameter a's of bound vortex panels (Pig, 7) as follows : 






^2 = 




- I 

0 




"x = + ”5 “^Sx 


Wfc = - a. 


I - a. - aj-S 
; 4 5 


= - / - a^OdS 




Similarly# 


and ^2 ” ^2^ ^3 ^ 



The velocity induced (Pig, 8 a) at any point C(x,y, 2) 
due to bound vortex is calculated using Biot-Savart law as 
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follows ; 
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[((y-y^) (2-22)-{y-y2) i-CCx-x^) (z-z^) ) 

j + ( (y-y2) - ^ “■^y*‘y2^ 

(z-z^) ) + ((x-x^) {Z-Z2) “ (-‘-^2^ (z-z^))^ + ( (x-x^) (y-y2) 

” (^-X2^ (y-Yl))^] {A-3) 


For a horse-shoe vortex 


From Fig, 8(b) V = + Z obA + 2- 
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APPENDIX - B 


Svaluation of inteQrals of inauced. veloci.t:y 
From Ec. (3,7) 
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where r = [(x-l + y^ + (z- = )^ ]^ 

+ a^^ 
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Subs ti tu ting 

(x- 1) = [y2+ (z- S 2 j^tan 6 = K tan 0 


and integrating 
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For calculating induced velocity due to i=l, =(^)= 0 is 
Eq, of line BC (Fig. 6b) (’)= 0 is the equation of AB in 

local coordinate system with origin at A. For i=2 origin is 
shifted to B and 0 0 are equations of CD and BC^ 

f: 

= - is the distance 
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Solving with the same substitution as (x-;^) = K tan 
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It is easy to apply simpson’s rule for evaluating these 
integrals rather than looking for closed form expressions. 
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Quadrilateral panels over the flat plate are not subdivided 
unlike free vortex panels over the wake, then summation from 
Eq. (A-5) is dropped and the equations ^^(^)= 0, 0 

becomes the Eq, of lines AS and CD respectively (Pig. 6b), 

Second integral of A-5 remains the same from 0 to 



59 


^PPSUDIX ~ C 

Boundary Conditions over bound vozt&z T^anels 

Ilie boundairy condition ecuations are satisfied at certain 
nodes of bound vortex panels to octain the undetermined 
coefficients of the local vorticity distribution over bound 
vortex panels. 


The no-penetration condition Ec. (2,6) is enfoi'ced at the 

* •f*]h 

average point of the bound vortex panels. At m control point 
Eg, (2.6) is of the form 
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Where ^ is the angle of attack and suffix m refers to 
the receiver panel and k to the sender panel, d # d 

til 

are the cosines of angle between normal at m control point 
and x,y,z axis respectively. 

For a flat plate rectangular wing if we are calculating 
the induced velocity due to bound vortex, (C-l) reduces to 

N 5 

£ E ^^j^m,k ^ “ sin a {C-2) 

ks=l j=l 
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hes to be evaluated 


I'/here d v;it:i suffix in the cosine of the angle betx'/een 
its two subscripts on sender panel and Uj, v^. , are the 
velocity components per unit . 

In addition we x-^rite equation of vorticity continuity at 

each cCOTnon node of bound vortex panels i,e, w„ and w„ should 

z X 

be continuous at each global node/ (Fig. 9 ) , 


w 

z 

11 


“S' 

d 

&2 

w 

X 

= ’'£'^5 X 

4- 

W| 

d 

IX 


For example 

w (1,4) = w (II/l) 

z z 

w (1/2) = w CII, 1) 

z z (C-3) 

w (III/ 3) = w (IV, 2) 
z z 

w (II/2) = (IV/4) 

z ^ 

and four equations fron continuity at common node (1,3), i.e. 

w (1,3) = w (II/2) 

3 Z 

w (III,4) = w (IV,1) 

z z (C-.4) 

w (1,3) = w (III,4) 

(11,2) » (IV, 1) 
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Continuity in x-direction is obtained by replacing z 
by X frcra (c-3) and (c-4) . 

Kutta condition, Eq, (2,10) is enforced at the nodes 
along the edges of separation, such as (III, 2), (III, 3), (IV, 2) 
& (P/, 3), Ref, Pig, 9 , where 

- w V =0 
X z z X 


cosp - sin p 


(c-s) 


where P is the angle with z-axis. For a flat plate 
rectangular wing this P is taken as 90'’, 

Kutta condition at these four nodes is 

w„ (III, 2 ) = 0 

w, (III, 3) = 0 

(IV, 2 ) = 0 

Z 

{IV, 3) = 0 

The last set of equations is obtained from symmetry 
condition for a symmetric flow along thie root chord, that 
vorticity in x-direction is zero, 

= 0 

= 

= Wfc for a flat plate rectangular wing 


{C- 6 ) 
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( 1 , 1 ) = 0 

( 1 , 2 ) =0 

(C-7) 

M (III, 1) = 0 
(III, 2) = 0 

The overdetermined set of stations obtained from (C—2) to 
iC—7) is solved by the method of least squares. 
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appendix - D 


Calculation of unkno-m« her: 

.coefficients of bound vortex oanela. 


-t fl'-’i -» * 


;2.S ‘tHZTj 


It IS difficult to solve the unknovms (vorticity distribution 
coefricients) of wake panels individually. If we satisfy continuity 
and symmetry condition over the v/ake panels it gives under- 
determined set of equations which cannot be solved. Therefore 
unknowns for free vortex panels are derived in terms of unknowns 
(a's) of bound vortex panels. This has been achieved by applying 
the condition that out flow of vorticity is equal to inflow 
of vorticity to free vortex panels (Pig. 10) . 


Outflovj of vorticity in z-direction = vorticity from panel 1 

^ vorticity from panel 5 

1 1 

= ” { (D-l) 


Inflow of vorticity = vorticity in z-direction to panel 9 


I 


1 


/ 

0 


(W^)g dZ 


(D-2) 


'"xh or 5 = + ''e 6x>l or 5 = or 5 

Equating (D-l) and (D-2) 


f (w^)g dz = / (Wg dz + / (Wg)^ dz 


0 


y* (a^+a^ ^ I ) d| =- / (a^ 1+ a^ + 5 ) d 


0 


^ (32 1+ + ®5 5)5^ 5 
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(a. 


“f 




t2 

*1 


3 2 '9 


). = 


/■?> fc / 

^ 2 ^ ’’1 


a , 4 
4 1 


+ a. 



^ 4- + 3. 

^ 1 4 1 




5 ^V<°-3> 


ii>iiic 0 flow leaves tlie TS at an angle of 90^# ^ 

a rectangular wing. 


t 

j> 

1 


for 


cosff ic^zionts of 4 ^ ^ front 


^^2^9 ~ ■" ^®2^5 


(D-4) 


■*■ “3 -^>9 = <®4 + =5 ^ >1 

At pt A, = 0 

and (a^)^ = -(a^)^ - 5 



+ ar 



(D-5) 


(D-6) 

{D-7) 


Similarly if we equate the outflow of vorticity from 
bound vor tex panels (1 and 5) to the inflow of vorticity to 
free vortex panel (9) in x— direction, 

we get, 

4 1 4 

/ <'"2^9 *= = S + } <'>2h ^ 

0 U u 

(w^)9 = Wj = (-a^t- )g 

or 5 = «? = + ^2?+ ^3* >1 or 5 


Prom(D-8) 



\fX 
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'4 “ ^ 5 ^ ^9 = i + ^ 2 ’ ^ + 


•3 = "1 


^1 

/ 


1 "** ^2 ^ "*" ^3 ^ ^ 5 ^ 


(D-9) 


Integrating 


(-a^?* ^ i~ ®4 = 1 “ ®5 




^1 ^2^ ^1 ■*■ ®3 2 ^5 (D-IO) 


Subtracting (D-4) from D-10) 


^^ 4 , ^5 


Ji, ) 

2 ^9 


®3 2 \ 


— (a* + a 


3“r ^5 


(D-ll) 


For a large aspect ratio wing at point A 

^ ^ = 0 
1 1 


(a.)— = — (a^ ) ^ •• (^^)c 
4 9 11 15 


(D-12) 


(35)^ = - (33)^^ - ^® 3^5 


(D-13) 


Outflow of vorticity from panel 9 eclated to inflow of 
vorticity to panel 10 can be written as follows 


/ (w^)9 dl= / (w^),o dl 
0 0 


(D-14) 


/ ^"^ 2^9 ~ S ^^ 2^10 

0 0 


(D-15) 
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{D-14) gives 


2 

1 


(a. . + a — i- ) 

11 2 13 2*^9 


|2 

(a^ >1+ 32" "1 + ^3 ”T ^IC 


(D~16) 


{^ 3)9 " ^^2^10 


I-, 


and (a^ +3-3 ^ \ = (a- + a_ — ~ ) _ 

2 ^9 1 3 2 ^10 


1 ‘ '^3 
and from D-15 




(-331 ^ ^ = (- a^^, r^- a^ 


(D--17) 

(D~1S) 


y2 

^ li % (D~19) 

5 2 ‘'10 


Prom D-17 and D-19 


5 


(a* + a_ ) sr (a + a — — ) 

4 5 2 '9 ' 4 ^ ^5 2 '’lO 


(D-20) 


Continuity condition at the points P and B (Pig. 10 ) in 
z~diroction gives 


(W,4 )q / 


P&B 


= <''^>10 / 


F&B 


(D-21) 


At point B (- - a^)^ = 

At point P (-a^ - = (a^ - a^-a^ 


Subtracting (D-22) from {D-23) 

(83 li - 5^)^ = (a 2 5^)^^ 

Fran D-17 and D-24 


10 


(D-22) 

(D-23) 

{D-24) 


(85)^ - (a5)^Q 


(D-25) 
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From D-24 and D-25 


^^4^ 1C 


(D— 26) 


Continuity in x-direction at 


points P and B rives 






iD-27) 


at point B (a^ + *]_)9 




2S ) 


at point F ^1^9 ~ ^-5 


1 " “2 ^ ^3 "I'lO 


Subtracting 


{D-29) 


(a^ $1 - a-3 ^ 2^9 ^^2 h " ^3 


(D-30) 


Prom D-17 and D-30 


( 33 )^ - (a3)3_o 


(D-31) 


Add D-27 to D-28 


we get (a^)^ = 


CD-32) 


Siniilarly vorticity distribution coefficients for panel 
10 and 17 are obtained by applying the continuity of vorticity 
in X and s-direction at points P and G (Pig. 10) 

continuity in z-direction 




(D-33) 


^ I ^10 ^17 / 


(D-34) 


pt G (a^ a^ - a^ ^ 


(D-35) 
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Subtracting 

"" ^^ 2^17 (D- 37 ) 

Continuity in x-dircction 

(Wy )in / = (Wy / ( 2 - 38 ) 

^ i.J-' -^y r2 


pt G 

pt F 


(o.^ + a_ Sa ^ 

1 2 1-^2 2.^10 


(a^ + ^2 ^ 1^10 == + ^3 


^*^1^ 17 
^ 1^17 


Subtracting 


^® 3^10 

(33)17 

Applying outflov 7 conditions 

we will get 



= (^1)17 



(33) 

= (23)17 


(D-39) 

(D-40) 


(D-41) 




z 


X 

Fig. 1 Arrangement of bound and free vortex panels. 



X 


No -Penetration Condition 
o Continuity of Vorticity Condition 
• Kutta Condition 
□ Symmetry Condition 

Kinematic and Dynamic Condition 

Fig. 2 Details of boundary conditions. 







Fig. 3 Quadrilateral and triangular vortex panels. 



Fig. 4 Spanwise circulation and shed vorticity 
distribution- 










>> z 



(D ® 

(2) (D 

n 

n 


X 

Fig. 9 Details of boundary conditions. 



Fig.lO Equating outflow of vorticity to inflow. 





















z/2b 


Fig. 12 Spanwise variation of pressure coefficient: AR=1.0, 
oC = 15deg , NX =5, NZ = 5 . 




Fig. 13 Convergence study fora rectangular wing: AR r 8,0 , oC a 5°, Xw/b=1.69 
NX =1 . NW =4. 













Hu Ho 


Fi^.17 Wake shape at different angle of attack for a rectangular wing : 
AR = 8.0, NX = 1, NZ =10, NW=7. 
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Z2tl)=2ih-NY) 

SO f U 79 - 

Y2h)ay3CI-l) 

ZZCDsZSiX^l) 

CONTINUE - 
DO 71 IJb 1,HW 
DO 71 JJal.NP 
Ittl EJ-1)^WP+J0 
JQstIi*«P 

irti.€x.so)Z2a>=7.|ci-SP) 
irci.sf •JO/zi ci^z^ct) 
irci,€T,J9}Z3CI>«Z2Cl) ^ _ 

YCI)«»|y 1 u5 + Y2{l5+Y'3C!i i<|. 

S C I J a«i h HI2 i I it 2:iy ll ^ ; 

0^5) 




U S| ' ' 

Mtr* ' 


n'CI),Y2ClJ#Yi€l)#YCl7 



'•f 1 €1 ) )»f 2f C Z3 C 1 3 -Zl CD) *^23 


ffYU3*»2+t«CI3**23 

S7 ' 

/f«(S3 


so 

Xlci3»t3CI*«P3 

Yta3*X3CJ-!«P3 



i«y 


111 

110 


42 

41 

28 


C 

C 

C 


c 

c 




/a U.i=ii3c;-Np) 

IMtis'll c u+sii»Tyt U/IKtl) 
iFt ro 111 
lF{l,F.i,l)X2C lj=xc 
XFt I.F:Q,1)//2CI}=0.U 

{^ ( C J ji^EQ ! 1 ?! A in! ? i ! ^ E, 1 1 ) X2 ( X ) XX u r-Mp+i ) 

1 F t C J,1 . KU , i ) , A«>), C 1 . i'lE , 1 ) ) X2 C r ) Xlf i c x-’jp+ u 

XECC JJ.Ea*l),A!MU,( l,NE,l)lii-«f(X)=ii3Cl-"^P+il 

Gi ta 110 

K2f l)rX3(I-l j 

5f2(£5»XUI-l5 

42(IJxZ?C I«U 

Xf llxCXl £I)+X2£I)+X3CI) )/3, 

XCU*fyUX> + X2tJj + £3Cl)l/J, 
idCD^CZlC £)+Z2£I) + ^3tl)l / 1. 

XM(I)=yi CIJ+CV3CI)-Xlf 
Z''JfI)xZlCI)<'CZ3{l)-ZltI J 

CONTINUE 
CO JTINUE 
CONTINUE 
STOP 
SNU 

SUBROUTINE CORi?( XM , YU , X J , T I , X2, Y2 , AN , N5f ) 

PROGRAM FOH FINDING THE POINTS OYER^WTNG SURfACE 

01 MgNSlUN XM(WY) .YMCMY) ^ XI (NY) , Yl (NX) , X2 ( NY ) , Ya C NY) , AN (N Y) 

i.SI.P£ 2 Cl),SLdPt 26 ) 

CpMMON/Ab/AUPHA.NX.NW 

OPfiNCuNIXxiJH^OEYICEx'OSK? .FHiEb^PUU.IM^) 
0PEW(UNIXs25,DEYICE»'05R,t,5 
READC38,») AR,SPAN,X|R,S^PUE 
SNPliExSNEEP AT U,R, IN DEGREES 
fPRxCR/CX 
PAIaS, 141592674 
XA»0,0 
YAw0,-O 

AREAB((SPAN)»*2)/ftR 

•t*5'«C0SDCSNp£4fi)t'2)5 



|C|04')£P)/€4,*NXO*CXA4‘C(«*I 
||||lc|D-KS)^C4,4NX3-£XA+f 


£X:A+C(N*»1,0 


l#f»Sl»PTIi)t(XC*»XA)/I4,f!«X3 + CCXC-XA34CR*l.))/^’liOAf €NX> 
f'l/FiOAtpf) 

YicofJsKpclHCXC*XA)/(4,4NX) + CCXC-'XA)*CN^t))/NX 
!a|McS3isEol(«5 + C^C"XA)l'3,/(4.»NX) + CCXC"XA)f CS*^1,3)/NX 



2 'i 

Af? 


C 

c 


c 

c 


'c' 

C 

c 


' , 1^1 

C’ 

C 


Cn^'TlMUE 
Cnf<T£M{J*2 
«»'■:£ UK*? 

KifMCl’TUN VKt.2f X i,X4,AM,yi, lf4,yrA,Zi,i'-‘A,2;M J 

ftsXfi-XJ 

a=:XM-X4 

C^yM^Yi 

l)=ifi'i-y4 

lf;aZN-ZJ 

F»Z’'i**Z4 

VEtiZa«CAA#0-ti»CJ/C (C=*:F-D^‘P:j»»2 + (A*F-B*E)’(‘#2+C^*0-«ff!)»*2) 
IJfCC (B-A)*A+CD-C)»C+tF-E)«'K)/SORT(A!)sA + C»n+E»E)-C CB-A)*B 
i + Cn-C>*U+CF-E)yp”)/SQKI’Ati*B+U*04F’‘'F) ) 

e«i) 

FU.'ICl'lQM VEliYCXi,X4,XN,y i,y4,yM,Z ?,Z4,Z^J) 

A»XN*XJ 

fcjs:XN-X4 

Csyw-y^ 

l)eyW.Y4 

e? 5 Zr** 7.3 

yECiyaCCS»E-A4Fj/CCC’»'F-D4S)»*2+CA4F-ti4E)442 + (A*U-B4C)4J2)) 
14tCC8-AJ*A+(f)-C)4C + (P’-e)4EJ/SQRlCA*A4C*C+E*E3-C CB-A)*B+ 
lC0-C)»r>+CF-E)*r3/SQRf (B48 + l>»r>4.F»F)) 

Ri:r»mN 

fir^CflOM TAMOC8E2A) 

rA!«i)BSlFn>CBETA)/C0S0C8ETA) 

KEfURM 

EMP 

FlNCflUM S|WP5M(A,8,H,.FJ 

»& ' 
ssCB-Rj/M 

S«sF'a)>r'Ut«'^W3;, 


'’ 1 ''': ''' ' ■ ■ ■ , ^ 

Mf YspCt/ C »*BtC*€ J ) 1* t i * O+A/SQRT C A*A ti»S+C4C J ) 
RErURK 
ENP 





c 

c 


c 

92 


1.01 


c 

c 


,i.r ’ ; 


c 

92 

101 


'.■ik 


ft ^ XN-Xb 
ii = yo-Vb 
C=ZM-Zb 
i / Lj ^ rzs - 
KFIilKfJ 

iiWD 


CB/CB>»'B+C’»'C))»Cl,OtA/SURl'f ft^'A + B*B + C’«=C) J 


FHNCrTQri CINTIRCXK) 

COMMLiri/Al/XX, Aft,BR,XP,E2,MM 
CDHM0N/A2/YYftrZZR 

«iKXT£;C2S/92)XX,YYR,ZZK,AA,BB,XP,E2,ClMTlR,XE,M/^ 

F0RHATC12X,9Fll .6,6X, 12) 
irCMM,E0,2)G0 10 tOl ^ ^ 

CMMatAA+BB*Xfe;)fZZ.R/C CXyH-XE)^»2+ZZR«*2) 

ZTE13-XP^XE/E2+XP 

Cl »iT| RaCMi^^: C Sr^ZTEl -SNZLKl ) 

KEIURO 

AAtS8»X«->BB»E2)»ZZS/( tYYR-XE)<'*2 + ZZR»*2) 
iolEi?aCclx-ltiE2)/SURT( CXX-Z£,e2)*J2 + (YYR-XE) J|2+ZZR»»2)) 

SBZTE2»C{,XX-XPJ/SQRTC CXX-Xp)*»2 + (XyR-XE)*#2 + ZZH#=i^2) > 

CINTlRoCHM#CSWZTE2-3MZLifc:2) 

RErURfI 

SMD 

niHcnon cintiLim ) 

cb«M0N/Al/XX,AA,0R;XP,e2,*!?(| 

eoMf«ari/Aj/fYC,,zzu 

§ fC « M . ea , 2 ) C»o I’D 101 
«?WcliAtOR»XE)*ZZ!4/C Cyytj-XE>**2+ZZli**2) 

Zf eia-XP»XE/E2fXP 

rqRHftTC16X,9Kll,6,&X,l2) 

C«Ms2cAA+BB*X€*§B*E2)»ZWiafU*3^E>**2<»'Zj;i.»»2) ■ 


' fk ' 





f 0 B f XP , i 2 » CIU 12 R , m , , Zf ei » S^SitiEl 

92 rolSllti | p 5 r 4 : rI 3 ) ' ' ■ 

a g j II & M / 

101 C»^MaCAA+l4B*rxe-8B'fe2)*CYYR-XE)/C(YYR-XE)*f2<l'Z;2,8*»2) 

ZLE 2 = XP ’»' Xe / E 2 



ton ^nn n n ^nn 


SNif>k.2=t txx-iil.fc:2 j/5«RT{, CXX-£t.F2)*’^2+ (YYR-Xh’)^=’t'2 + 27,K»*2) ) 
S''.'2'i'l':? = ( CXX-Xf’J/SQRTC CXX-XPJ*^24 cyyR-Xt.;)+^t2 + ZZH^»2) ) 
CII^'T'2R=CM^^»CS'!?.TE:;2-SM2I:jK2) 

^JRITeC23,RiUX, YYH,ZZH, AA,8B,Xp,£2,CIW'r2R,XK,Cf^W ZLf!:2^Sf’JZIje2 
3 ,S^^ZTii2>M 

3 FORWATC J iP’y.4,133 

e;wo 


2 

01 


93 


22 

23 


FUWCTION CiNT2uCXE3 
CO)4MUM/Al/XX,AA^Bn,XP,£2,M'1 

CQman/AS/i^u,zzh 

JlFCHM,Ea,2>GU ru 1U1 

CH»«=-CAAtBH»XE)*cy YL-XEJ/( CYYr<-XE)3t4=24-ZZb»f23 
ZfeJa-XP3tXE/K2tXP 

S'^Zl.ElaCXX/SQKf (XX»»2 + CYYt-XE)+*2 + ZZlb*»2) ) „ ^ ^ 

SWZTEi = C CXX-ZTeU/SQRTC CXX-ZTEJ 3 3'*2 + CyYL-XE)*»2 + ZZL’('f2)‘» 
Cl!«T2l4=CMM*(SNZTeJ -SNZBEl) 

<<RlTE(23/92JXX,yYL,ZZL.-AA,BB,Xp,E2,CIMT2li,XE,CMM,ZTei,SHZl4El 

1 ,SMZTE1,MM 

F0R^AT(2X,HFy.4,T2) 

HEfUKN 

Cf<(i«aCAA + BB#XE-BBtF2J3t( YYL-XEJ/C f YyLi*-XE)=l'»2 + Z2rLf *2) 
ZLE2aXP*XE/E2 , « 

Sf^ZEEZst tXX-ZLE2)/S«KT( (XX-Zt,£2)*»2 + CyYCi-XE)»*2 + ZZE*»23 3 
SWZT£2«C UX-XP3/S0RT( CXX-Xp3»*2+cyyo-XE)**2 + ZZli*»23 > 
eiMT2LsCf<M»CSNZTE2-SMZbE2) ^ „ „ „ „ 

«fRifSC23,93)XX,yyri,ZZL,AA,Btt,Xp,E2,C|NT2B,XE,CMX,Z|,E2f SNZEE2 

rgl‘^ATC2i,i3ry,'i,i23 

aiyuRn* 

l»l«30flNe G£MIRycXl,X2,X3,XM,fl,Y2#Y3,yp<,l«r) 

9I^ES??10R XlCWf >,yi(Nr3,X2CWr),y2CJir3 ,X3(Rr),y3CNF3 ^XMCNP 
13,YM€Nf3 

Ct1M#^Ofj/A6/A0PHA,.1X,f««« 

C0WM0r»|/A7/wY 
Y«=9*5 
WZsN Y+l 

3EfA«0,0 , , ^ 

ZZs9O,/F0OAT€f«X>' ^ ■ , 

tAilastnO'CifirAjXCO^BCftW') ' ' 

:gSi ■ ' ' ■ ", 

'IliilSJllfi,-? ■ ■ . .' 

fMf , 







{rp.ST,iSx3c«*^|BL09<»-o,6#a-i5 

XMC«|aCXlCM}+XaCH3tX3(H3i/3* 
Y'^Ci^3aaiCM3 + Y2CM3+y3tRJ3/3* 



111 


C 

c 

i4 

H*> 

95 


99 




20 


K2(;'''lj3X2C'^)+l}.i25 
X3f-nsXi(.'0 ^0.125 

co.rrxN'ie: ■ 
co.rrxNU!': 

KKlMiar? 

&Nl) 

SOtJROUTlNE EINVIA^N, AX,MJ 

ai4Ei4f>XnN AIN,H) , AICRr'O 
■42=riti 

00 i4 ijlsl,-''! 

00 14 l4J = M2,H 

ACtii.r,j)3‘),o 

on R5 K = lri‘l 

I2SKI-N 

ACK, 

DO 99 00=1, N 
J2bL J+ i 
psacoOfOO) 

DO 95 1=1, M 
ACOJ,I)aACLJ,IJ/P 
00 99 LK=1,N 
00 99 Ot»J2,M 
IFCLK-LJ ) 20, 99, 2U 

A ( U , 01 ) =A C fjK , Cl 1 1 - A ( liO , Oi ) t A ( OK , L 0 ) 

COfiTlNUE 
i>0 100 1 = 1, H 
00 100 J=M2,M 

;':|f(I,03)=Atl,J> 

I'.'RBifyRN 
' UMy 

5lKBr«Crj5?^Wi!:?^?"n5|0),,«DExp50,2),D|n50) 

EilUtJstElicE c{RO»,5R5»);(ICnbU«,jCDfcu'l),(««AX,T,SK4P) 

l>EieR« = 1,0 
do 20 J=1,N 
IP£VOT(J)=0,a 

do 550 1*1,*’ ' ' ' , ' 

A>^AX=O*0 , , ■ 



40; ^ 


ACt'C5l!#5S|01fcSi«Af 
irCM) 260,260/210 
do 250 0=1, 

BCiRoSJ'0}e»!icOOyM,l.) 



25y dCiCUf.U«/4j=SWAP 
ib'J i-iUPKC I, i > = IRUw 

tf-'iTex.t l,2) = T.C0Lnn 

f> I n i‘ s A c I ca L a M , I'CHL'i >» ) 
D'’’ix)^=pTyrn’ 
ki lCQr.U!S^cnLimj=i .n 
do 35U L=i,M 

i50 ACiCafjU"<,»45sA(ICOr,lJfJ!,L)/pIVOf 
IP-(M) 3Rd,iRO,i6D 
36\) do 370 L= 1 ,M 

370 3CIC0LUM,I.J=BC XCULUM , t.) / PI VO f 

380 do S30 Li!=A,r'i 

il*-(0l-icur.u«) 400,5*^0,400 
400 r-ACbi ,icnouM) 

ACPI, £COLUM)=:0,0 
do 450 f.sl.fj 

450 ACH ,b)=AlLl,0)-A(lC(JLUM,b) 41‘ 
lF£d) 550,550,460 
460 d0'5y0 L=i,5 

500 d(Ll ,0IxBCC,l,lr)-BCXC0LUM,0J4r 

550 CndTIOUF: 

do 710 1 = 1 , N 
OxW+l"! 


0Ei’p:i?w=f>eTfe:RM4aTC0) 

IFC rfJDKKdi, 1) -INDEX CL.,2J) 630,710,630 

J8a»I=IN0EKlt.,U 

3Ca6U«JSlN»>EXCb,2) 

do 7&§ 

SWAP«ACK,a«Ofe») 

ACK,3gOw5*5ACK,JCOLUM} 

A|K,jC!Ot.UM)aS*>AP 
CpfI«UE 
CWflNUE 
Jo it KsltN 

irtlPltfOTfKl.NE.l) GO TO U 
COWTINUE 


RCTHRN 


12 <iRiTE(5,y91) 

991 FORWATCloX, 'fUTKIX IS SISGUIiAR 
74‘i RETURN 


END 



S J IRnilTX UfiS') A® C A , B, 4 , .'i A , , I ft , 13 , I DGT, WK ARKA , ) 

nrc---- — »-S---- — -r.IHRARY 


i»3AGF: 

PARA'IiTriRS 




P 

,', -tiA’^aaAGg 


H 

I'i A 
iurt 

lA 

IB 

tnsT 

*|KAREA 

IRR 


enuTiNgs 


) 


LEAST SafJARgS SDLUrinN OF nvEHDE’^ERs^IMRD 
SYSrEB OF LiHEAH EQUATIONS 
call LLS'^ARC MA,NP, TA,IB, Tonr, i^KAREAr ter 

T.JF CDEFFICIEI^P MATRIX 6f TUE ' EQfl ATI OM 

AX = B, WHERE A IS M X f-.A WITH M GREATER I'HA'^ 

'iR Ey:jAL r'T ihput a is replaced by 

THE P3EUD3-1 ’‘’^EHSE dF A 

Matrix of the right ra'^d side if the eouattom 


AX = Bf where B is 'I X .\B, THE MA K ?ia 
SDLjTiOfj X aVERwlRlTES B, 

NUBBER OF ROWS IM ft ftNO B, 

■'lU'^BER OF CDLuMf^S TR MATRIX A 

WUB3ER OS' roLUMr^S Ji\ matrix b 
ROW niHENSTOM OF A IM THE CALLING 
PROGRaB, 

ROW DIMEMSTOM of b in the calling 
program, 

the elements of a ARff: assumed to be correct 

TO IDGT SlG;nFICANT DIGITS. IDGT IS AN 
TNPJT parameter, 

WORK AREA OF OEHENSION GREATER THAN OR EOUAL 
TO NA(NA+4). 

ERROR PARAMfiTgR 

terminal error s 128 + N 

Nat INDICATES THAT INPUT A IS A EERO 
MATRIX, 

SINGLE 

LPSaOR^LSVALR^HERTST, VSORTM 


mm iimmm mmmm #r jpt wyi* pi 



* IT p f a a A " . j i A i\ ^ s ' 'i 

[.•niiHAr - iRi-’A*'' RiJ-! rr.ilf-T SilM'-lAftV 


rMPnRF'iNJl'; iTor 
rnutine see tf^e 
// ... // n^V’ of* 


a cfininiete s*5erif iCrat ion ot 
t'AH Fni<rpaiy [(lOrarv Manual 
3 moif'ifientat Ion dppen<^ent. 


the use of t*^is 
Terms mar*<e'l 


A,_^Purpose 


//f t>i A AP*/ caiC'iiatos the aoproxifnate inverse of a real matrix 
hv Croat's metnod. 


B, Specification 


C 

c 


S'HJUHTI.v 


K //s'ul AA(?'// 

ITAiTU) 

ca, h 

//real// AUA 


, TBNIT 
,H), iJMi 


( A , I A , , 

iFATf 

T( lU'JlT^N) , 


OfjTr, TUNIT, WKSPCS, 
WKSpCECrtl 


C, Parameters 

Kssassassss^c 

h • //real// array of 0 TME?viSinN (lA#p? wnere p,gE#M, 

eefore entry, k mast contain the elements ot the real 

VacceS'Sf ai exit, it contains the Ctout factorisation with 
tn'r anitddiagonai of n unaerstooa, 

lA - T^jfBOSP. 

On entry, lA must sneclfy the first dimension of array A as 
declared in the calling tsiiblnrogram*. 
lA.GE.fl, , ' . , , 

unchanged on exit# . 

U • IS'fRSSK, , ! '' , , 

Oh': entry# ;S, 'hi 'fftfi^trix a* 

imitfpy Where o , SE . N . : 
the inverse of 




..•o.ii-ir» . ...I.* «p*cl»y the first aieenslon of arrasf 

,,.J?S|’ 4 S<!tir'e 9 *Jfi thS calUns CsebJPtoara*. 

unchanged on exit, 

WKSPCE - //real// array of DlwewsiQN at least (n)* 


Used as sorting space. 



iP’Ari. 


i.-IfKURK 


8Ptorf> «ntrv, 1^’nlli Mns*" be asst-jne'^ n v^lue, for users not 
l^iTitllisr vHb this parameter (ffescribei tn Chaoter POJ in 
tne 'iiatn iaG P’fiR'j'RAH Lihrarv 'Manual) tne recofnff'ended value 
is 

unless the routine detects an error (see next section), 
l'5^Arij contains <) on exit. 

u. Error Tnilcators and ‘■kernings 
Errors aetectea nv tne routine;- 


l^’AfU = 1 

j’oe natrix A is sinouiar or altnost singular. Possibly due 
to rounain) errors. 


£*^'D OF P^IAAF FHrtTRAfj SOMMAPY - MARK 10 
FEBRUARY 1‘^Ri 

♦ f"EH0 OF F')l 



